<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html>
 <head>
  <title>cbdsqr.f</title>
 <meta name="generator" content="emacs 21.3.1; htmlfontify 0.20">
<style type="text/css"><!-- 
body { background: rgb(255, 255, 255);  color: rgb(0, 0, 0);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.default   { background: rgb(255, 255, 255);  color: rgb(0, 0, 0);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.default a { background: rgb(255, 255, 255);  color: rgb(0, 0, 0);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: underline; }
span.string   { color: rgb(188, 143, 143);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.string a { color: rgb(188, 143, 143);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: underline; }
span.comment   { color: rgb(178, 34, 34);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.comment a { color: rgb(178, 34, 34);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: underline; }
 --></style>

 </head>
  <body>

<pre>
      SUBROUTINE <a name="CBDSQR.1"></a><a href="cbdsqr.f.html#CBDSQR.1">CBDSQR</a>( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U,
     $                   LDU, C, LDC, RWORK, INFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  -- LAPACK routine (version 3.1) --
</span><span class="comment">*</span><span class="comment">     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
</span><span class="comment">*</span><span class="comment">     November 2006
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     .. Scalar Arguments ..
</span>      CHARACTER          UPLO
      INTEGER            INFO, LDC, LDU, LDVT, N, NCC, NCVT, NRU
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Array Arguments ..
</span>      REAL               D( * ), E( * ), RWORK( * )
      COMPLEX            C( LDC, * ), U( LDU, * ), VT( LDVT, * )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Purpose
</span><span class="comment">*</span><span class="comment">  =======
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  <a name="CBDSQR.20"></a><a href="cbdsqr.f.html#CBDSQR.1">CBDSQR</a> computes the singular values and, optionally, the right and/or
</span><span class="comment">*</span><span class="comment">  left singular vectors from the singular value decomposition (SVD) of
</span><span class="comment">*</span><span class="comment">  a real N-by-N (upper or lower) bidiagonal matrix B using the implicit
</span><span class="comment">*</span><span class="comment">  zero-shift QR algorithm.  The SVD of B has the form
</span><span class="comment">*</span><span class="comment">  
</span><span class="comment">*</span><span class="comment">     B = Q * S * P**H
</span><span class="comment">*</span><span class="comment">  
</span><span class="comment">*</span><span class="comment">  where S is the diagonal matrix of singular values, Q is an orthogonal
</span><span class="comment">*</span><span class="comment">  matrix of left singular vectors, and P is an orthogonal matrix of
</span><span class="comment">*</span><span class="comment">  right singular vectors.  If left singular vectors are requested, this
</span><span class="comment">*</span><span class="comment">  subroutine actually returns U*Q instead of Q, and, if right singular
</span><span class="comment">*</span><span class="comment">  vectors are requested, this subroutine returns P**H*VT instead of
</span><span class="comment">*</span><span class="comment">  P**H, for given complex input matrices U and VT.  When U and VT are
</span><span class="comment">*</span><span class="comment">  the unitary matrices that reduce a general matrix A to bidiagonal
</span><span class="comment">*</span><span class="comment">  form: A = U*B*VT, as computed by <a name="CGEBRD.34"></a><a href="cgebrd.f.html#CGEBRD.1">CGEBRD</a>, then
</span><span class="comment">*</span><span class="comment">  
</span><span class="comment">*</span><span class="comment">     A = (U*Q) * S * (P**H*VT)
</span><span class="comment">*</span><span class="comment">  
</span><span class="comment">*</span><span class="comment">  is the SVD of A.  Optionally, the subroutine may also compute Q**H*C
</span><span class="comment">*</span><span class="comment">  for a given complex input matrix C.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  See &quot;Computing  Small Singular Values of Bidiagonal Matrices With
</span><span class="comment">*</span><span class="comment">  Guaranteed High Relative Accuracy,&quot; by J. Demmel and W. Kahan,
</span><span class="comment">*</span><span class="comment">  LAPACK Working Note #3 (or SIAM J. Sci. Statist. Comput. vol. 11,
</span><span class="comment">*</span><span class="comment">  no. 5, pp. 873-912, Sept 1990) and
</span><span class="comment">*</span><span class="comment">  &quot;Accurate singular values and differential qd algorithms,&quot; by
</span><span class="comment">*</span><span class="comment">  B. Parlett and V. Fernando, Technical Report CPAM-554, Mathematics
</span><span class="comment">*</span><span class="comment">  Department, University of California at Berkeley, July 1992
</span><span class="comment">*</span><span class="comment">  for a detailed description of the algorithm.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Arguments
</span><span class="comment">*</span><span class="comment">  =========
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  UPLO    (input) CHARACTER*1
</span><span class="comment">*</span><span class="comment">          = 'U':  B is upper bidiagonal;
</span><span class="comment">*</span><span class="comment">          = 'L':  B is lower bidiagonal.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  N       (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The order of the matrix B.  N &gt;= 0.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  NCVT    (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The number of columns of the matrix VT. NCVT &gt;= 0.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  NRU     (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The number of rows of the matrix U. NRU &gt;= 0.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  NCC     (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The number of columns of the matrix C. NCC &gt;= 0.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  D       (input/output) REAL array, dimension (N)
</span><span class="comment">*</span><span class="comment">          On entry, the n diagonal elements of the bidiagonal matrix B.
</span><span class="comment">*</span><span class="comment">          On exit, if INFO=0, the singular values of B in decreasing
</span><span class="comment">*</span><span class="comment">          order.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  E       (input/output) REAL array, dimension (N-1)
</span><span class="comment">*</span><span class="comment">          On entry, the N-1 offdiagonal elements of the bidiagonal
</span><span class="comment">*</span><span class="comment">          matrix B.
</span><span class="comment">*</span><span class="comment">          On exit, if INFO = 0, E is destroyed; if INFO &gt; 0, D and E
</span><span class="comment">*</span><span class="comment">          will contain the diagonal and superdiagonal elements of a
</span><span class="comment">*</span><span class="comment">          bidiagonal matrix orthogonally equivalent to the one given
</span><span class="comment">*</span><span class="comment">          as input.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  VT      (input/output) COMPLEX array, dimension (LDVT, NCVT)
</span><span class="comment">*</span><span class="comment">          On entry, an N-by-NCVT matrix VT.
</span><span class="comment">*</span><span class="comment">          On exit, VT is overwritten by P**H * VT.
</span><span class="comment">*</span><span class="comment">          Not referenced if NCVT = 0.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  LDVT    (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The leading dimension of the array VT.
</span><span class="comment">*</span><span class="comment">          LDVT &gt;= max(1,N) if NCVT &gt; 0; LDVT &gt;= 1 if NCVT = 0.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  U       (input/output) COMPLEX array, dimension (LDU, N)
</span><span class="comment">*</span><span class="comment">          On entry, an NRU-by-N matrix U.
</span><span class="comment">*</span><span class="comment">          On exit, U is overwritten by U * Q.
</span><span class="comment">*</span><span class="comment">          Not referenced if NRU = 0.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  LDU     (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The leading dimension of the array U.  LDU &gt;= max(1,NRU).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  C       (input/output) COMPLEX array, dimension (LDC, NCC)
</span><span class="comment">*</span><span class="comment">          On entry, an N-by-NCC matrix C.
</span><span class="comment">*</span><span class="comment">          On exit, C is overwritten by Q**H * C.
</span><span class="comment">*</span><span class="comment">          Not referenced if NCC = 0.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  LDC     (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The leading dimension of the array C.
</span><span class="comment">*</span><span class="comment">          LDC &gt;= max(1,N) if NCC &gt; 0; LDC &gt;=1 if NCC = 0.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  RWORK   (workspace) REAL array, dimension (2*N) 
</span><span class="comment">*</span><span class="comment">          if NCVT = NRU = NCC = 0, (max(1, 4*N-4)) otherwise
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  INFO    (output) INTEGER
</span><span class="comment">*</span><span class="comment">          = 0:  successful exit
</span><span class="comment">*</span><span class="comment">          &lt; 0:  If INFO = -i, the i-th argument had an illegal value
</span><span class="comment">*</span><span class="comment">          &gt; 0:  the algorithm did not converge; D and E contain the
</span><span class="comment">*</span><span class="comment">                elements of a bidiagonal matrix which is orthogonally
</span><span class="comment">*</span><span class="comment">                similar to the input matrix B;  if INFO = i, i
</span><span class="comment">*</span><span class="comment">                elements of E have not converged to zero.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Internal Parameters
</span><span class="comment">*</span><span class="comment">  ===================
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  TOLMUL  REAL, default = max(10,min(100,EPS**(-1/8)))
</span><span class="comment">*</span><span class="comment">          TOLMUL controls the convergence criterion of the QR loop.
</span><span class="comment">*</span><span class="comment">          If it is positive, TOLMUL*EPS is the desired relative
</span><span class="comment">*</span><span class="comment">             precision in the computed singular values.
</span><span class="comment">*</span><span class="comment">          If it is negative, abs(TOLMUL*EPS*sigma_max) is the
</span><span class="comment">*</span><span class="comment">             desired absolute accuracy in the computed singular
</span><span class="comment">*</span><span class="comment">             values (corresponds to relative accuracy
</span><span class="comment">*</span><span class="comment">             abs(TOLMUL*EPS) in the largest singular value.
</span><span class="comment">*</span><span class="comment">          abs(TOLMUL) should be between 1 and 1/EPS, and preferably
</span><span class="comment">*</span><span class="comment">             between 10 (for fast convergence) and .1/EPS
</span><span class="comment">*</span><span class="comment">             (for there to be some accuracy in the results).
</span><span class="comment">*</span><span class="comment">          Default is to lose at either one eighth or 2 of the
</span><span class="comment">*</span><span class="comment">             available decimal digits in each computed singular value
</span><span class="comment">*</span><span class="comment">             (whichever is smaller).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  MAXITR  INTEGER, default = 6
</span><span class="comment">*</span><span class="comment">          MAXITR controls the maximum number of passes of the
</span><span class="comment">*</span><span class="comment">          algorithm through its inner loop. The algorithms stops
</span><span class="comment">*</span><span class="comment">          (and so fails to converge) if the number of passes
</span><span class="comment">*</span><span class="comment">          through the inner loop exceeds MAXITR*N**2.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  =====================================================================
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     .. Parameters ..
</span>      REAL               ZERO
      PARAMETER          ( ZERO = 0.0E0 )
      REAL               ONE
      PARAMETER          ( ONE = 1.0E0 )
      REAL               NEGONE
      PARAMETER          ( NEGONE = -1.0E0 )
      REAL               HNDRTH
      PARAMETER          ( HNDRTH = 0.01E0 )
      REAL               TEN
      PARAMETER          ( TEN = 10.0E0 )
      REAL               HNDRD
      PARAMETER          ( HNDRD = 100.0E0 )
      REAL               MEIGTH
      PARAMETER          ( MEIGTH = -0.125E0 )
      INTEGER            MAXITR
      PARAMETER          ( MAXITR = 6 )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Scalars ..
</span>      LOGICAL            LOWER, ROTATE
      INTEGER            I, IDIR, ISUB, ITER, J, LL, LLL, M, MAXIT, NM1,
     $                   NM12, NM13, OLDLL, OLDM
      REAL               ABSE, ABSS, COSL, COSR, CS, EPS, F, G, H, MU,
     $                   OLDCS, OLDSN, R, SHIFT, SIGMN, SIGMX, SINL,
     $                   SINR, SLL, SMAX, SMIN, SMINL, SMINOA,
     $                   SN, THRESH, TOL, TOLMUL, UNFL
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Functions ..
</span>      LOGICAL            <a name="LSAME.173"></a><a href="lsame.f.html#LSAME.1">LSAME</a>
      REAL               <a name="SLAMCH.174"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>
      EXTERNAL           <a name="LSAME.175"></a><a href="lsame.f.html#LSAME.1">LSAME</a>, <a name="SLAMCH.175"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Subroutines ..
</span>      EXTERNAL           <a name="CLASR.178"></a><a href="clasr.f.html#CLASR.1">CLASR</a>, <a name="CSROT.178"></a><a href="csrot.f.html#CSROT.1">CSROT</a>, CSSCAL, CSWAP, <a name="SLARTG.178"></a><a href="slartg.f.html#SLARTG.1">SLARTG</a>, <a name="SLAS2.178"></a><a href="slas2.f.html#SLAS2.1">SLAS2</a>,
     $                   <a name="SLASQ1.179"></a><a href="slasq1.f.html#SLASQ1.1">SLASQ1</a>, <a name="SLASV2.179"></a><a href="slasv2.f.html#SLASV2.1">SLASV2</a>, <a name="XERBLA.179"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Intrinsic Functions ..
</span>      INTRINSIC          ABS, MAX, MIN, REAL, SIGN, SQRT
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Executable Statements ..
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Test the input parameters.
</span><span class="comment">*</span><span class="comment">
</span>      INFO = 0
      LOWER = <a name="LSAME.189"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( UPLO, <span class="string">'L'</span> )
      IF( .NOT.<a name="LSAME.190"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( UPLO, <span class="string">'U'</span> ) .AND. .NOT.LOWER ) THEN
         INFO = -1
      ELSE IF( N.LT.0 ) THEN
         INFO = -2
      ELSE IF( NCVT.LT.0 ) THEN
         INFO = -3
      ELSE IF( NRU.LT.0 ) THEN
         INFO = -4
      ELSE IF( NCC.LT.0 ) THEN
         INFO = -5
      ELSE IF( ( NCVT.EQ.0 .AND. LDVT.LT.1 ) .OR.
     $         ( NCVT.GT.0 .AND. LDVT.LT.MAX( 1, N ) ) ) THEN
         INFO = -9
      ELSE IF( LDU.LT.MAX( 1, NRU ) ) THEN
         INFO = -11
      ELSE IF( ( NCC.EQ.0 .AND. LDC.LT.1 ) .OR.
     $         ( NCC.GT.0 .AND. LDC.LT.MAX( 1, N ) ) ) THEN
         INFO = -13
      END IF
      IF( INFO.NE.0 ) THEN
         CALL <a name="XERBLA.210"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="CBDSQR.210"></a><a href="cbdsqr.f.html#CBDSQR.1">CBDSQR</a>'</span>, -INFO )
         RETURN
      END IF
      IF( N.EQ.0 )
     $   RETURN
      IF( N.EQ.1 )
     $   GO TO 160
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     ROTATE is true if any singular vectors desired, false otherwise
</span><span class="comment">*</span><span class="comment">
</span>      ROTATE = ( NCVT.GT.0 ) .OR. ( NRU.GT.0 ) .OR. ( NCC.GT.0 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     If no singular vectors desired, use qd algorithm
</span><span class="comment">*</span><span class="comment">
</span>      IF( .NOT.ROTATE ) THEN
         CALL <a name="SLASQ1.225"></a><a href="slasq1.f.html#SLASQ1.1">SLASQ1</a>( N, D, E, RWORK, INFO )
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span>      NM1 = N - 1
      NM12 = NM1 + NM1
      NM13 = NM12 + NM1
      IDIR = 0
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Get machine constants
</span><span class="comment">*</span><span class="comment">
</span>      EPS = <a name="SLAMCH.236"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>( <span class="string">'Epsilon'</span> )
      UNFL = <a name="SLAMCH.237"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>( <span class="string">'Safe minimum'</span> )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     If matrix lower bidiagonal, rotate to be upper bidiagonal
</span><span class="comment">*</span><span class="comment">     by applying Givens rotations on the left
</span><span class="comment">*</span><span class="comment">
</span>      IF( LOWER ) THEN
         DO 10 I = 1, N - 1
            CALL <a name="SLARTG.244"></a><a href="slartg.f.html#SLARTG.1">SLARTG</a>( D( I ), E( I ), CS, SN, R )
            D( I ) = R
            E( I ) = SN*D( I+1 )
            D( I+1 ) = CS*D( I+1 )
            RWORK( I ) = CS
            RWORK( NM1+I ) = SN
   10    CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Update singular vectors if desired
</span><span class="comment">*</span><span class="comment">
</span>         IF( NRU.GT.0 )
     $      CALL <a name="CLASR.255"></a><a href="clasr.f.html#CLASR.1">CLASR</a>( <span class="string">'R'</span>, <span class="string">'V'</span>, <span class="string">'F'</span>, NRU, N, RWORK( 1 ), RWORK( N ),
     $                  U, LDU )
         IF( NCC.GT.0 )
     $      CALL <a name="CLASR.258"></a><a href="clasr.f.html#CLASR.1">CLASR</a>( <span class="string">'L'</span>, <span class="string">'V'</span>, <span class="string">'F'</span>, N, NCC, RWORK( 1 ), RWORK( N ),
     $                  C, LDC )
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Compute singular values to relative accuracy TOL
</span><span class="comment">*</span><span class="comment">     (By setting TOL to be negative, algorithm will compute
</span><span class="comment">*</span><span class="comment">     singular values to absolute accuracy ABS(TOL)*norm(input matrix))
</span><span class="comment">*</span><span class="comment">
</span>      TOLMUL = MAX( TEN, MIN( HNDRD, EPS**MEIGTH ) )
      TOL = TOLMUL*EPS
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Compute approximate maximum, minimum singular values
</span><span class="comment">*</span><span class="comment">
</span>      SMAX = ZERO
      DO 20 I = 1, N
         SMAX = MAX( SMAX, ABS( D( I ) ) )
   20 CONTINUE
      DO 30 I = 1, N - 1
         SMAX = MAX( SMAX, ABS( E( I ) ) )
   30 CONTINUE
      SMINL = ZERO
      IF( TOL.GE.ZERO ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Relative accuracy desired
</span><span class="comment">*</span><span class="comment">
</span>         SMINOA = ABS( D( 1 ) )
         IF( SMINOA.EQ.ZERO )
     $      GO TO 50
         MU = SMINOA
         DO 40 I = 2, N
            MU = ABS( D( I ) )*( MU / ( MU+ABS( E( I-1 ) ) ) )
            SMINOA = MIN( SMINOA, MU )
            IF( SMINOA.EQ.ZERO )
     $         GO TO 50
   40    CONTINUE
   50    CONTINUE
         SMINOA = SMINOA / SQRT( REAL( N ) )
         THRESH = MAX( TOL*SMINOA, MAXITR*N*N*UNFL )
      ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Absolute accuracy desired
</span><span class="comment">*</span><span class="comment">
</span>         THRESH = MAX( ABS( TOL )*SMAX, MAXITR*N*N*UNFL )
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Prepare for main iteration loop for the singular values
</span><span class="comment">*</span><span class="comment">     (MAXIT is the maximum number of passes through the inner
</span><span class="comment">*</span><span class="comment">     loop permitted before nonconvergence signalled.)
</span><span class="comment">*</span><span class="comment">
</span>      MAXIT = MAXITR*N*N
      ITER = 0
      OLDLL = -1
      OLDM = -1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     M points to last element of unconverged part of matrix
</span><span class="comment">*</span><span class="comment">
</span>      M = N
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Begin main iteration loop
</span><span class="comment">*</span><span class="comment">
</span>   60 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Check for convergence or exceeding iteration count
</span><span class="comment">*</span><span class="comment">
</span>      IF( M.LE.1 )
     $   GO TO 160
      IF( ITER.GT.MAXIT )
     $   GO TO 200
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Find diagonal block of matrix to work on
</span><span class="comment">*</span><span class="comment">
</span>      IF( TOL.LT.ZERO .AND. ABS( D( M ) ).LE.THRESH )
     $   D( M ) = ZERO
      SMAX = ABS( D( M ) )
      SMIN = SMAX
      DO 70 LLL = 1, M - 1
         LL = M - LLL
         ABSS = ABS( D( LL ) )
         ABSE = ABS( E( LL ) )
         IF( TOL.LT.ZERO .AND. ABSS.LE.THRESH )
     $      D( LL ) = ZERO
         IF( ABSE.LE.THRESH )
     $      GO TO 80
         SMIN = MIN( SMIN, ABSS )
         SMAX = MAX( SMAX, ABSS, ABSE )
   70 CONTINUE
      LL = 0
      GO TO 90
   80 CONTINUE
      E( LL ) = ZERO
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Matrix splits since E(LL) = 0
</span><span class="comment">*</span><span class="comment">
</span>      IF( LL.EQ.M-1 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Convergence of bottom singular value, return to top of loop
</span><span class="comment">*</span><span class="comment">
</span>         M = M - 1
         GO TO 60
      END IF
   90 CONTINUE
      LL = LL + 1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     E(LL) through E(M-1) are nonzero, E(LL-1) is zero
</span><span class="comment">*</span><span class="comment">
</span>      IF( LL.EQ.M-1 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        2 by 2 block, handle separately
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="SLASV2.367"></a><a href="slasv2.f.html#SLASV2.1">SLASV2</a>( D( M-1 ), E( M-1 ), D( M ), SIGMN, SIGMX, SINR,
     $                COSR, SINL, COSL )
         D( M-1 ) = SIGMX
         E( M-1 ) = ZERO
         D( M ) = SIGMN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Compute singular vectors, if desired
</span><span class="comment">*</span><span class="comment">
</span>         IF( NCVT.GT.0 )
     $      CALL <a name="CSROT.376"></a><a href="csrot.f.html#CSROT.1">CSROT</a>( NCVT, VT( M-1, 1 ), LDVT, VT( M, 1 ), LDVT,
     $                  COSR, SINR )
         IF( NRU.GT.0 )
     $      CALL <a name="CSROT.379"></a><a href="csrot.f.html#CSROT.1">CSROT</a>( NRU, U( 1, M-1 ), 1, U( 1, M ), 1, COSL, SINL )
         IF( NCC.GT.0 )
     $      CALL <a name="CSROT.381"></a><a href="csrot.f.html#CSROT.1">CSROT</a>( NCC, C( M-1, 1 ), LDC, C( M, 1 ), LDC, COSL,
     $                  SINL )
         M = M - 2
         GO TO 60
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     If working on new submatrix, choose shift direction
</span><span class="comment">*</span><span class="comment">     (from larger end diagonal element towards smaller)
</span><span class="comment">*</span><span class="comment">
</span>      IF( LL.GT.OLDM .OR. M.LT.OLDLL ) THEN
         IF( ABS( D( LL ) ).GE.ABS( D( M ) ) ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Chase bulge from top (big end) to bottom (small end)
</span><span class="comment">*</span><span class="comment">
</span>            IDIR = 1
         ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Chase bulge from bottom (big end) to top (small end)
</span><span class="comment">*</span><span class="comment">
</span>            IDIR = 2
         END IF
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Apply convergence tests
</span><span class="comment">*</span><span class="comment">
</span>      IF( IDIR.EQ.1 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Run convergence test in forward direction
</span><span class="comment">*</span><span class="comment">        First apply standard test to bottom of matrix
</span><span class="comment">*</span><span class="comment">
</span>         IF( ABS( E( M-1 ) ).LE.ABS( TOL )*ABS( D( M ) ) .OR.
     $       ( TOL.LT.ZERO .AND. ABS( E( M-1 ) ).LE.THRESH ) ) THEN
            E( M-1 ) = ZERO
            GO TO 60
         END IF
<span class="comment">*</span><span class="comment">
</span>         IF( TOL.GE.ZERO ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           If relative accuracy desired,
</span><span class="comment">*</span><span class="comment">           apply convergence criterion forward
</span><span class="comment">*</span><span class="comment">
</span>            MU = ABS( D( LL ) )
            SMINL = MU
            DO 100 LLL = LL, M - 1
               IF( ABS( E( LLL ) ).LE.TOL*MU ) THEN
                  E( LLL ) = ZERO
                  GO TO 60
               END IF
               MU = ABS( D( LLL+1 ) )*( MU / ( MU+ABS( E( LLL ) ) ) )
               SMINL = MIN( SMINL, MU )
  100       CONTINUE
         END IF
<span class="comment">*</span><span class="comment">
</span>      ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Run convergence test in backward direction
</span><span class="comment">*</span><span class="comment">        First apply standard test to top of matrix
</span><span class="comment">*</span><span class="comment">
</span>         IF( ABS( E( LL ) ).LE.ABS( TOL )*ABS( D( LL ) ) .OR.
     $       ( TOL.LT.ZERO .AND. ABS( E( LL ) ).LE.THRESH ) ) THEN
            E( LL ) = ZERO
            GO TO 60
         END IF
<span class="comment">*</span><span class="comment">
</span>         IF( TOL.GE.ZERO ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           If relative accuracy desired,
</span><span class="comment">*</span><span class="comment">           apply convergence criterion backward
</span><span class="comment">*</span><span class="comment">
</span>            MU = ABS( D( M ) )
            SMINL = MU
            DO 110 LLL = M - 1, LL, -1
               IF( ABS( E( LLL ) ).LE.TOL*MU ) THEN
                  E( LLL ) = ZERO
                  GO TO 60
               END IF
               MU = ABS( D( LLL ) )*( MU / ( MU+ABS( E( LLL ) ) ) )
               SMINL = MIN( SMINL, MU )
  110       CONTINUE
         END IF
      END IF
      OLDLL = LL
      OLDM = M
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Compute shift.  First, test if shifting would ruin relative
</span><span class="comment">*</span><span class="comment">     accuracy, and if so set the shift to zero.
</span><span class="comment">*</span><span class="comment">
</span>      IF( TOL.GE.ZERO .AND. N*TOL*( SMINL / SMAX ).LE.
     $    MAX( EPS, HNDRTH*TOL ) ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Use a zero shift to avoid loss of relative accuracy
</span><span class="comment">*</span><span class="comment">
</span>         SHIFT = ZERO
      ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Compute the shift from 2-by-2 block at end of matrix
</span><span class="comment">*</span><span class="comment">
</span>         IF( IDIR.EQ.1 ) THEN
            SLL = ABS( D( LL ) )
            CALL <a name="SLAS2.480"></a><a href="slas2.f.html#SLAS2.1">SLAS2</a>( D( M-1 ), E( M-1 ), D( M ), SHIFT, R )
         ELSE
            SLL = ABS( D( M ) )
            CALL <a name="SLAS2.483"></a><a href="slas2.f.html#SLAS2.1">SLAS2</a>( D( LL ), E( LL ), D( LL+1 ), SHIFT, R )
         END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Test if shift negligible, and if so set to zero
</span><span class="comment">*</span><span class="comment">
</span>         IF( SLL.GT.ZERO ) THEN
            IF( ( SHIFT / SLL )**2.LT.EPS )
     $         SHIFT = ZERO
         END IF
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Increment iteration count
</span><span class="comment">*</span><span class="comment">
</span>      ITER = ITER + M - LL
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     If SHIFT = 0, do simplified QR iteration
</span><span class="comment">*</span><span class="comment">
</span>      IF( SHIFT.EQ.ZERO ) THEN
         IF( IDIR.EQ.1 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Chase bulge from top to bottom
</span><span class="comment">*</span><span class="comment">           Save cosines and sines for later singular vector updates
</span><span class="comment">*</span><span class="comment">
</span>            CS = ONE
            OLDCS = ONE
            DO 120 I = LL, M - 1
               CALL <a name="SLARTG.509"></a><a href="slartg.f.html#SLARTG.1">SLARTG</a>( D( I )*CS, E( I ), CS, SN, R )
               IF( I.GT.LL )
     $            E( I-1 ) = OLDSN*R
               CALL <a name="SLARTG.512"></a><a href="slartg.f.html#SLARTG.1">SLARTG</a>( OLDCS*R, D( I+1 )*SN, OLDCS, OLDSN, D( I ) )
               RWORK( I-LL+1 ) = CS
               RWORK( I-LL+1+NM1 ) = SN
               RWORK( I-LL+1+NM12 ) = OLDCS
               RWORK( I-LL+1+NM13 ) = OLDSN
  120       CONTINUE
            H = D( M )*CS
            D( M ) = H*OLDCS
            E( M-1 ) = H*OLDSN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Update singular vectors
</span><span class="comment">*</span><span class="comment">
</span>            IF( NCVT.GT.0 )
     $         CALL <a name="CLASR.525"></a><a href="clasr.f.html#CLASR.1">CLASR</a>( <span class="string">'L'</span>, <span class="string">'V'</span>, <span class="string">'F'</span>, M-LL+1, NCVT, RWORK( 1 ),
     $                     RWORK( N ), VT( LL, 1 ), LDVT )
            IF( NRU.GT.0 )
     $         CALL <a name="CLASR.528"></a><a href="clasr.f.html#CLASR.1">CLASR</a>( <span class="string">'R'</span>, <span class="string">'V'</span>, <span class="string">'F'</span>, NRU, M-LL+1, RWORK( NM12+1 ),
     $                     RWORK( NM13+1 ), U( 1, LL ), LDU )
            IF( NCC.GT.0 )
     $         CALL <a name="CLASR.531"></a><a href="clasr.f.html#CLASR.1">CLASR</a>( <span class="string">'L'</span>, <span class="string">'V'</span>, <span class="string">'F'</span>, M-LL+1, NCC, RWORK( NM12+1 ),
     $                     RWORK( NM13+1 ), C( LL, 1 ), LDC )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Test convergence
</span><span class="comment">*</span><span class="comment">
</span>            IF( ABS( E( M-1 ) ).LE.THRESH )
     $         E( M-1 ) = ZERO
<span class="comment">*</span><span class="comment">
</span>         ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Chase bulge from bottom to top
</span><span class="comment">*</span><span class="comment">           Save cosines and sines for later singular vector updates
</span><span class="comment">*</span><span class="comment">
</span>            CS = ONE
            OLDCS = ONE
            DO 130 I = M, LL + 1, -1
               CALL <a name="SLARTG.547"></a><a href="slartg.f.html#SLARTG.1">SLARTG</a>( D( I )*CS, E( I-1 ), CS, SN, R )
               IF( I.LT.M )
     $            E( I ) = OLDSN*R
               CALL <a name="SLARTG.550"></a><a href="slartg.f.html#SLARTG.1">SLARTG</a>( OLDCS*R, D( I-1 )*SN, OLDCS, OLDSN, D( I ) )
               RWORK( I-LL ) = CS
               RWORK( I-LL+NM1 ) = -SN
               RWORK( I-LL+NM12 ) = OLDCS
               RWORK( I-LL+NM13 ) = -OLDSN
  130       CONTINUE
            H = D( LL )*CS
            D( LL ) = H*OLDCS
            E( LL ) = H*OLDSN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Update singular vectors
</span><span class="comment">*</span><span class="comment">
</span>            IF( NCVT.GT.0 )
     $         CALL <a name="CLASR.563"></a><a href="clasr.f.html#CLASR.1">CLASR</a>( <span class="string">'L'</span>, <span class="string">'V'</span>, <span class="string">'B'</span>, M-LL+1, NCVT, RWORK( NM12+1 ),
     $                     RWORK( NM13+1 ), VT( LL, 1 ), LDVT )
            IF( NRU.GT.0 )
     $         CALL <a name="CLASR.566"></a><a href="clasr.f.html#CLASR.1">CLASR</a>( <span class="string">'R'</span>, <span class="string">'V'</span>, <span class="string">'B'</span>, NRU, M-LL+1, RWORK( 1 ),
     $                     RWORK( N ), U( 1, LL ), LDU )
            IF( NCC.GT.0 )
     $         CALL <a name="CLASR.569"></a><a href="clasr.f.html#CLASR.1">CLASR</a>( <span class="string">'L'</span>, <span class="string">'V'</span>, <span class="string">'B'</span>, M-LL+1, NCC, RWORK( 1 ),
     $                     RWORK( N ), C( LL, 1 ), LDC )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Test convergence
</span><span class="comment">*</span><span class="comment">
</span>            IF( ABS( E( LL ) ).LE.THRESH )
     $         E( LL ) = ZERO
         END IF
      ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Use nonzero shift
</span><span class="comment">*</span><span class="comment">
</span>         IF( IDIR.EQ.1 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Chase bulge from top to bottom
</span><span class="comment">*</span><span class="comment">           Save cosines and sines for later singular vector updates
</span><span class="comment">*</span><span class="comment">
</span>            F = ( ABS( D( LL ) )-SHIFT )*
     $          ( SIGN( ONE, D( LL ) )+SHIFT / D( LL ) )
            G = E( LL )
            DO 140 I = LL, M - 1
               CALL <a name="SLARTG.590"></a><a href="slartg.f.html#SLARTG.1">SLARTG</a>( F, G, COSR, SINR, R )
               IF( I.GT.LL )
     $            E( I-1 ) = R
               F = COSR*D( I ) + SINR*E( I )
               E( I ) = COSR*E( I ) - SINR*D( I )
               G = SINR*D( I+1 )
               D( I+1 ) = COSR*D( I+1 )
               CALL <a name="SLARTG.597"></a><a href="slartg.f.html#SLARTG.1">SLARTG</a>( F, G, COSL, SINL, R )
               D( I ) = R
               F = COSL*E( I ) + SINL*D( I+1 )
               D( I+1 ) = COSL*D( I+1 ) - SINL*E( I )
               IF( I.LT.M-1 ) THEN
                  G = SINL*E( I+1 )
                  E( I+1 ) = COSL*E( I+1 )
               END IF
               RWORK( I-LL+1 ) = COSR
               RWORK( I-LL+1+NM1 ) = SINR
               RWORK( I-LL+1+NM12 ) = COSL
               RWORK( I-LL+1+NM13 ) = SINL
  140       CONTINUE
            E( M-1 ) = F
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Update singular vectors
</span><span class="comment">*</span><span class="comment">
</span>            IF( NCVT.GT.0 )
     $         CALL <a name="CLASR.615"></a><a href="clasr.f.html#CLASR.1">CLASR</a>( <span class="string">'L'</span>, <span class="string">'V'</span>, <span class="string">'F'</span>, M-LL+1, NCVT, RWORK( 1 ),
     $                     RWORK( N ), VT( LL, 1 ), LDVT )
            IF( NRU.GT.0 )
     $         CALL <a name="CLASR.618"></a><a href="clasr.f.html#CLASR.1">CLASR</a>( <span class="string">'R'</span>, <span class="string">'V'</span>, <span class="string">'F'</span>, NRU, M-LL+1, RWORK( NM12+1 ),
     $                     RWORK( NM13+1 ), U( 1, LL ), LDU )
            IF( NCC.GT.0 )
     $         CALL <a name="CLASR.621"></a><a href="clasr.f.html#CLASR.1">CLASR</a>( <span class="string">'L'</span>, <span class="string">'V'</span>, <span class="string">'F'</span>, M-LL+1, NCC, RWORK( NM12+1 ),
     $                     RWORK( NM13+1 ), C( LL, 1 ), LDC )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Test convergence
</span><span class="comment">*</span><span class="comment">
</span>            IF( ABS( E( M-1 ) ).LE.THRESH )
     $         E( M-1 ) = ZERO
<span class="comment">*</span><span class="comment">
</span>         ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Chase bulge from bottom to top
</span><span class="comment">*</span><span class="comment">           Save cosines and sines for later singular vector updates
</span><span class="comment">*</span><span class="comment">
</span>            F = ( ABS( D( M ) )-SHIFT )*( SIGN( ONE, D( M ) )+SHIFT /
     $          D( M ) )
            G = E( M-1 )
            DO 150 I = M, LL + 1, -1
               CALL <a name="SLARTG.638"></a><a href="slartg.f.html#SLARTG.1">SLARTG</a>( F, G, COSR, SINR, R )
               IF( I.LT.M )
     $            E( I ) = R
               F = COSR*D( I ) + SINR*E( I-1 )
               E( I-1 ) = COSR*E( I-1 ) - SINR*D( I )
               G = SINR*D( I-1 )
               D( I-1 ) = COSR*D( I-1 )
               CALL <a name="SLARTG.645"></a><a href="slartg.f.html#SLARTG.1">SLARTG</a>( F, G, COSL, SINL, R )
               D( I ) = R
               F = COSL*E( I-1 ) + SINL*D( I-1 )
               D( I-1 ) = COSL*D( I-1 ) - SINL*E( I-1 )
               IF( I.GT.LL+1 ) THEN
                  G = SINL*E( I-2 )
                  E( I-2 ) = COSL*E( I-2 )
               END IF
               RWORK( I-LL ) = COSR
               RWORK( I-LL+NM1 ) = -SINR
               RWORK( I-LL+NM12 ) = COSL
               RWORK( I-LL+NM13 ) = -SINL
  150       CONTINUE
            E( LL ) = F
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Test convergence
</span><span class="comment">*</span><span class="comment">
</span>            IF( ABS( E( LL ) ).LE.THRESH )
     $         E( LL ) = ZERO
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Update singular vectors if desired
</span><span class="comment">*</span><span class="comment">
</span>            IF( NCVT.GT.0 )
     $         CALL <a name="CLASR.668"></a><a href="clasr.f.html#CLASR.1">CLASR</a>( <span class="string">'L'</span>, <span class="string">'V'</span>, <span class="string">'B'</span>, M-LL+1, NCVT, RWORK( NM12+1 ),
     $                     RWORK( NM13+1 ), VT( LL, 1 ), LDVT )
            IF( NRU.GT.0 )
     $         CALL <a name="CLASR.671"></a><a href="clasr.f.html#CLASR.1">CLASR</a>( <span class="string">'R'</span>, <span class="string">'V'</span>, <span class="string">'B'</span>, NRU, M-LL+1, RWORK( 1 ),
     $                     RWORK( N ), U( 1, LL ), LDU )
            IF( NCC.GT.0 )
     $         CALL <a name="CLASR.674"></a><a href="clasr.f.html#CLASR.1">CLASR</a>( <span class="string">'L'</span>, <span class="string">'V'</span>, <span class="string">'B'</span>, M-LL+1, NCC, RWORK( 1 ),
     $                     RWORK( N ), C( LL, 1 ), LDC )
         END IF
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     QR iteration finished, go back and check convergence
</span><span class="comment">*</span><span class="comment">
</span>      GO TO 60
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     All singular values converged, so make them positive
</span><span class="comment">*</span><span class="comment">
</span>  160 CONTINUE
      DO 170 I = 1, N
         IF( D( I ).LT.ZERO ) THEN
            D( I ) = -D( I )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Change sign of singular vectors, if desired
</span><span class="comment">*</span><span class="comment">
</span>            IF( NCVT.GT.0 )
     $         CALL CSSCAL( NCVT, NEGONE, VT( I, 1 ), LDVT )
         END IF
  170 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Sort the singular values into decreasing order (insertion sort on
</span><span class="comment">*</span><span class="comment">     singular values, but only one transposition per singular vector)
</span><span class="comment">*</span><span class="comment">
</span>      DO 190 I = 1, N - 1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Scan for smallest D(I)
</span><span class="comment">*</span><span class="comment">
</span>         ISUB = 1
         SMIN = D( 1 )
         DO 180 J = 2, N + 1 - I
            IF( D( J ).LE.SMIN ) THEN
               ISUB = J
               SMIN = D( J )
            END IF
  180    CONTINUE
         IF( ISUB.NE.N+1-I ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Swap singular values and vectors
</span><span class="comment">*</span><span class="comment">
</span>            D( ISUB ) = D( N+1-I )
            D( N+1-I ) = SMIN
            IF( NCVT.GT.0 )
     $         CALL CSWAP( NCVT, VT( ISUB, 1 ), LDVT, VT( N+1-I, 1 ),
     $                     LDVT )
            IF( NRU.GT.0 )
     $         CALL CSWAP( NRU, U( 1, ISUB ), 1, U( 1, N+1-I ), 1 )
            IF( NCC.GT.0 )
     $         CALL CSWAP( NCC, C( ISUB, 1 ), LDC, C( N+1-I, 1 ), LDC )
         END IF
  190 CONTINUE
      GO TO 220
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Maximum number of iterations exceeded, failure to converge
</span><span class="comment">*</span><span class="comment">
</span>  200 CONTINUE
      INFO = 0
      DO 210 I = 1, N - 1
         IF( E( I ).NE.ZERO )
     $      INFO = INFO + 1
  210 CONTINUE
  220 CONTINUE
      RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     End of <a name="CBDSQR.740"></a><a href="cbdsqr.f.html#CBDSQR.1">CBDSQR</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>
